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We present a tutorial for Bayesian inference of time-evolving coupled systems in the presence of noise. It 
includes the necessary theoretical description and the algorithms for its implementation. For general program- 
ming purposes, a pseudocode description is given. Examples based on coupled phase and limit-cycle oscillators 
illustrate the most important features. Codes written in MatLab for the method and the examples accompany 
the tutorial. 



I. INTRODUCTION 

The Bayesian methods for dynamical inference make for 
important signal processing techniques. Their applications 
include physics, biology and climate. Theoretical studies of 
such methods have already been presented extensively JU-ID]. 
Of particular interest is the method [j]} |H published recently, 
which can identify time-evolving coupled dynamics and is 
able to follow the time-variability of coupling functions. Here 
we discuss how this method can be implemented, including 
the algorithms, programming and applications. The problems 
and phenomena that are treated by the method are already rel- 
atively complex, hence in order to make the presentation more 
understandable we provide a different presentation style with 
examples that are simple and elementary. 

The tutorial is organized as follows. First, in section HI1 we 
present the theoretical terms needed for the implementation 
of Bayesian inference. The algorithms and the programming 
description are given in section[nH Application of the method 
on two examples is shown in section [IV] The first example 
uses coupled phase oscillators to present the basics of the in- 
ference of time-evolving phase dynamics in presence of noise. 
The second example uses coupled limit-cycle oscillators, pre- 
senting also the reconstruction and detection of synchroniza- 
tion and coupling functions. Some important generalizations 
of the framework are presented in section [V] Finally, short 
concluding remarks are given in section [VTl 



II. INFERENCE IMPLEMENTATION 

Given that I vectors of A^-samples time-series X = {<fri, n = 
t>i{t n )} (t n = nh) are provided, the main task for the method 



is to infer the unknown model parameters and the noise diffu- 



sion matrix M 



{4° 



Eij, }. The model for the interacting 



phase dynamics to be inferred (e.g. for I = 2 oscillators) is 
given by the following stochastic differential equations: 
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$2,0 = 



$u:(0i,0 2 ) + 6W, 
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1, Cq = u>i are the respective fre- 



where $ 1-0 

quencies, and the rest of and c[ are the K most im- 
portant Fourier components which serve as base functions for 
the inference. The noise is assumed to be white, Gaussian, 
and parameterized by a noise diffusion matrix (E 2 x 2): 

Umj(r)) = S(t-T)E ij . 

The inference exploits Bayes' theorem, assumes Euler mid- 
point discretisation [4)i tn = {(j)i, n+ i - 4>i, n )/h and 4>* l n = 
{4>i,n + 0;.ri+i)/2], and uses the properties of the stochas- 
tic integral to evaluate the likelihood function (see |Q]-[5l] for 
details). Assuming that the parameters c have a multivari- 
ate normal distribution, with mean c, and a covariance matrix 
S = S _1 , the stationary distribution is calculated recursively 
using only four equations: for the parameters c, for the noise 
matrix E, for the concentration matrix H and for temporary 
equation r, which will be defined bellow. These four equa- 
tions are applied on every time-series window, and they are 
central for the implementation of the inference. For easier im- 
plementation, and as an alternative to the theoretical presenta- 
tion HHH], in the following we present the matrix notation of 
these equations. 

For sake of clarification, ignoring the stochastic term, we 
can write EqQ]at every point in time t n in matrix form: 
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The big matrix in the center can be seen as composed by K 
blocks of diagonal /-sized matrices. We can define the latter 
matrix multiplication as: 

4>{t n ) = P(>„) x c. (2) 

From Eq. and the discussion in JH, the noise matrix E can 
be expressed as: 

h N T 

E= (^n)-P(tn)c) (J>(t n ) - P(t„)c) . (3) 

n=l 

Next, the equation for the concentration matrix H (i.e. the in- 
verse of the covariance matrix S) can be calculated as: 

N 

3 = H pnor + h ]T P(t n ) T E- 1 P(t„). (4) 

n=l 

Note that, while more compact, we do not suggest that the ma- 
trix notation should be implemented as it is because the matrix 
P is mostly filled with zeros and such implementation might 
produce very inefficient code. Similarly, one can express the 
temporary matrix variable r: 

N , N 

r = S prior c + h^2 P(t„) T E- 1 0(t n ) - - v(t„), (5) 

where each element on the lines of the vector v(t n ) = 
d<S>l d& ' s ' m practice, the elements of each column of 
P(i„) derived with respect of the l-th variable, where I is the 
corresponding line of the element of P(i„). 

Finally, after the results from the other calculations the pa- 
rameters vector c is evaluated as: 

c = E- 1 r. (6) 

By evaluating the four equations (fSJl-© using the readout 
time series X = {<fri, n = <t>i{tn)} one can effectively calcu- 
late the multivariate probability Mx (c|, c, H) which explicitly 
defines the probability density of each parameter set of the dy- 
namical system. Following the matrix notation and the respec- 
tive dimensions, one can implement the method consistently. 
For example, if we are given two phase time-series (I = 2) and 
we use six base functions in total (K = 6), then the respective 
dimensions of the matrices are: Ef X ; = E2 X 2, ?Kxi = ?6xi, 

BiKxK = 26x6, CixK = Ci X Q, PlxK = P2X6 and VA'xl = 
V 6 xl- 

The inference method needs to follow the time-evolution 
of the parameter set c while separating the effects from the 
noise. In order to achieve this we modify the propagation pro- 
cedure between the covariance of the current posterior £ post 

and the next prior E",^ 1 IB- The definite matrix E^g is intro- 
duced and shows how much each parameter diffuses normally. 
Thus, the next prior probability of the parameters is the convo- 
lution of two current normal multivariate distributions, S post 
and S diff : = S° st + ££ ff . We consider S diff to be 

such that there is no change of correlation between parame- 
ters (pij = 5ij) and that each standard deviation Oi is a known 



fraction of the relevant standard deviation from the posterior 
covariance (or parameters) Oi = p w (<r" ost )i, where p w is a 
constant parameter. In practice this means that Sjjff has zero 
values everywhere, except for the diagonal values, which are 
a fraction of the diagonal values of the posterior S post . 

III. ALGORITHMS AND PROGRAMMING 

In this section we discuss the algorithmic and program- 
ming details for the implementation. First we present arguably 
the most complicated part; the algorithm for the dynamical 
Bayesian inference applied on a single window of readout 
data. The algorithm employing recursion using the Eqs. (0- 
© can be summarized with the following steps: 

i) the algorithm starts from a c pr i or and E pr j 01 -, 

ii) noise matrix E new is calculated using Eq. (01, 

iii) H new is calculated using Eq. 

iv) r is calculated using Eq. ©, 

v) c new is calculated using Eq. (O, 

vi) then again to point ii) using c ne „ as c. 

The stopping rule is when a 'convergence' is reached i.e. 
when further running the algorithm would not modify c 
and H any further. For example, we used the condition: 
X) i c oid — c new) 2 /c^ew < £ where e is some very small con- 
stant. Since the problem is parabolic, this convergence is very 
fast - typically a few cycles J27j]. The prior distribution is as- 
sumed to be a noninformative "flat" acting as an initial limit 
of an infinitely large normal distribution, by setting H pr j or = 
and c prior = 0. 

For general programming purposes, in the follow- 
ing we outline an informal pseudo-code descrip- 
tion of the main algorithms and sub-algorithms. 
The comments are presented in grey. First we de- 
scribe the algorithm for the Bayesian inference: 



Algorithm 1 : Bayesian inference 

\\calculate temporary variables beforehand 
FORk=l:N 

- calculate P 

- calculate v 
FOR i=l:k 

FORj=l:k 

hPTP(i,j)=hPTP(i,j)+P(i mod l,i div Z)*P(j mod l,j div I) 
ENDFOR 
FOR g=l:Z 

hPtx(i,g)=hPtx(i,g)+P(imoJ I, i div l)*(t> g {k) 
ENDFOR 
ENDFOR 
hV2=hV2+v 
ENDFOR 

hPtx=h*hPtx; hPTP=h*hPTP; hV2=0.5*h*hV2 
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FOR lp= 1 :MaxLoops \\main recursive loop 

- calculate E 

- calculate c pt 

IF SUM((c pr - c pt ) 2 /c2 t )<SmallError 

RETURN 
ENDIF 

ENDFOR 



The standard functions mod and div are used only to position 
the calculation on the right indexes. Since indexes can have 
different notations in various programming languages (i.e. 
starting from zero i — 0, 1, ... or one i = 1,2,.. .), one 
should take care when applying mod and div in this context. 
The sub-algorithms 'calculate P' and 'calculate v' depend 
on the particular base functions and their partial derivatives. 
They only involve a simply evaluation of the base functions 
in respect to <j>*(k), and will not be discussed in detail. The 
other main calculations are performed within 'calculate E' 
and 'calculate c pt \ which are discussed in detail bellow. 

Algorithm 2: calculate E \\use of Eq. (0 
FOR i=l:N 
- calculate P 

E = E+(</> r P*c)*(transpose of (<&-P*c)) \\where <j) 2X N 
ENDFOR 
E = h/N*E 



Finally, the algorithm for calculation of the parameters is 
expressed as: 



Algorithm 3: calculate c 
invE = inverse of E 



The three described algorithms are applied to a single 
window of data. The time-series are separated into se- 
quential blocks of data and the algorithms are applied to 
each of them. The core of Bayesian inference is that the 
evaluation of the next block of data depends and uses the 
results of the evaluation from the previous block. The 
process of information propagation, between the posterior 
and the next prior distribution, can be adjusted to allow the 
time-variability of the parameters to be followed. We used 
propagation depending on the concentration matrix H pt or 
the parameters vector c pt . In the following we describe the 
propagation with respect to the concentration matrix H pt : 



Algorithm 4: Propagation 
invH" = inverse of H pt ; 
invDiff n =0 
FOR i=l:K 

invDiff"(i,i)=^*invH n (i,i) 
ENDFOR 

S" +1 =inverse of (invDiff"+invH") 



Once the inference is performed, one can use the inferred 
parameters to detect certain dynamical and phenomenological 
characteristics of the interacting systems. For example, calcu- 
lating the norm of the parameters which are inferred with the 
relevant base functions one can detect the coupling strength 
and directionality between the oscillators floT-llOll. Similarly, 
one can reconstruct the form of the coupling function Jll fill — 
[Till . To evaluate if the systems are synchronized, one can test 
if systems with such inferred parameters are undergoing syn- 
chronization lfl5l - il7ll . This can be done, for example, with 
maps representation and the modified Newton root-finding 
method (see I2] for details). 



\\calculate H, Eq. © 




\\calculate c, Eq. © 

c=(inverse of H pt )*r Figure 1: The instantaneous phases generated by the model of two 

coupled phase oscillators ([7}. (a) presents the phase from <f>i and (b) 
from 4>2- For simpler presentation, the phases are 'wrapped' in 2-7T. 



4 



IV. EXAMPLES 

This section outlines two examples for phase reconstruction 
of coupled oscillators with time-varying dynamics and subject 
to noise. The first example illustrates the basics of the infer- 
ence on a simple phase oscillators model, while the second 
example involves limit-cycle oscillators and presents the de- 
tection of several characteristics and relationships that can be 
detected from the inferred parameters. MatLab codes for the 
examples are provided at [28]. 



A. Coupled phase oscillators 

In order to present in a transparent way the basics of the 
inference technique we first consider two coupled phase os- 
cillators lfl8tl subject to noise: 



01 = wi(i) +aisin(0i) + a 3 (t) sin(0 2 ) + 6(0 

02 = w 2 + a 2 sin(0i) + a 4 sin(0 2 ) + 6(0- 



(7) 



Each oscillator is described by the frequency parameter ui, 
u>2, the parameters for the self dynamics ai, a 4 and the cou- 
pling parameters a 2 , a.3 for the direct influence from the other 
oscillator. Two parameters are set to be periodically time- 
varying, the frequency u>i(t) = 2 — 0.5 sin(27r0.00151f) and 
the coupling parameter 0,3(0 = 0.8 — 0.3 sin(27r0.0012£). 
The noises are set to be white Gaussian and uncorrected be- 
tween themselves. The rest of the parameters are uj 2 = 4.53, 
01 = 0.8, a 2 = 0, a 4 = 0.6, E u = 0.03 and E 22 = 0.01. 
Fig. Q] shows the nature and the noise perturbations of the 
instantaneous phases on which the Bayesian inference is ap- 
plied. 

The choice to use phase oscillators is very convenient for 
the inference because one needs to reconstruct the same phase 
dynamics. Therefore, the phase model is known beforehand 
and the deterministic terms of the rhs of coupled system (|7]i 
are the actual base functions to be used in order to infer the 
six parameters (cJi, W2, a\, a 2 , 0,3 and a 4 ). The results from 
the inference of one block of data are presented in Table [IV] 
One can notice that the method reconstructs the intrinsic pa- 
rameters successfully and with high precision. Additionally, 
and this is unique for this method, the level and the corre- 
lations of the noise are inferred very precisely. Fig. [2] shows 
the time-variations of all the parameters inferred from sequen- 
tial windows of length w — 40 s, and propagation constant 
p w = 0.2. It is easy to notice that the parameters and their 
time-variability are precisely inferred. 



B. Coupled limit-cycle oscillators 

The second example involves a system of two coupled 
limit-cycle oscillators, which can serve as a model for a 
number of oscillating processes in nature, including cardio- 
respiratory, electrochemical, mechanical, biological [U 
dHH^]. The model consists of two interacting Poincare os- 
cillators subject to noise: 
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Figure 2: Time-evolution of the inferred parameters from model Q. 
(a) and (b) present the two frequencies, (c) and (f) the self dynamics 
parameters, and (d) and (e) present the coupling parameters. The 
original intrinsic values of the parameters are presented with black 
color, underneath the inferred parameters. 



*i = - {\l x l + v\- ^) x i + £i(£2 -xi) + 6(0 

m = -(\jA + vi- 1)2/1 + +£1(2/2 -yi) + 6(0 



x 2 = - {\] x 2 + V2 - l) x 2 ~u 2 y 2 + e 2 (t)(xi -x 2 ) +6(0 
2/2 = -(\J x l + yl - 1)2/2 + ^2x2 + £2(0(2/1 - 2/2) + 6(0. 



(8) 



where periodic time-variability is introduced in the frequency 
of the first oscillator a>i(i) = 1 - 0.4sin(27r0.002t) and 
the coupling parameter from the first to the second oscillator 
e 2 (0 = 0.2-0.1 sin(27r0.00m). The noises are white Gaus- 
sian again, with no correlations between them. The rest of the 
parameters are u 2 = 4.91, £i = 0.05, En = E 22 = 0.007 
and E33 = E44 = 0.004. The systems are simulated in state 
space, and the corresponding signals and phase portrait of the 
first system are given in Fig. [3](a) and (b). 

The phases are estimated as </>j = arctan(j/i/a;i) (arctan 
being a four-quadrant function^ from the state signals. Alter- 
natively, one can use Hilbert B20I1 or synchrosqueezed trans- 
form lull . The choice of base functions for the inference 
needs to be determined so that the phase dynamics can be 
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Parameters 




W2 | 


ai \ 


a 2 | 


as | 


ai | 


£n| 


£12 | 


E 21 \ 


E22 


Intrinsic values 


2.032 


4.53 


0.8 





1.013 


0.6 


0.03 








0.01 


Inferred mean values 


2.026 


4.537 


0.803 


-0.014 


1.054 


0.596 


0.029 


-0.000 


-0.000 


0.010 



Table I: Results from the inference of numerically simulated system l(7J. The first row describes the physical meaning of the parameters, the 
last two rows show the values of the intrinsic parameters and their inferred mean values, respectively. The results are presented for one window 
of data around t — 1980 s. 
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Figure 3: Inference of time-varying parameters from model of two 
coupled limit-cycle oscillators l[8}. Phase portrait showing the noisy 
state of the first oscillator (a). The time-series of the two oscillators 
xi(t) and X2(t) (b) . The inferred time-evolution of the frequency 
parameter (c), and the net coupling from the first to the second oscil- 
lator (d). 



reconstructed effectively. Because of the oscillatory nature 
and the periodic solutions, we decompose the phase dynamics 
in Fourier series. Hence, the Fourier series up to some or- 
der serve as base functions for the Bayesian inference. Care 
must be taken to ensure that the functions are not very lin- 
early dependent, such as sin(x) and sin(— x), which can lead 
to imprecise and wrong separation of parameters within the 
inference. In the following, we used only one side of the ex- 
pansion, e.g. for sm(ncj>i + mfo) we used the components 
n = 1, . . . , K instead of n = —K, . . . ,K. The results for the 
reconstruction of the time-varying parameters are presented 
in Fig. [3(c) and (d). The periodic sine variations are evident 
both in the frequency and the coupling strength. Note that the 
coupling amplitude is evaluated as a norm of all the relevant 
inferred parameters that describe this influence. 

One can use the inferred parameters not only for evaluation 
of characteristics of individual oscillators, but also to deter- 
mine if the coupled system undergoes some qualitative tran- 
sitions. An obvious example of the latter is the detection of 
synchronization. For this reason we modify the parameters 
of model ® with wi(t) = 1 - 0.4sin(27r0.0022t), constant 
coupling parameter £2 = 0.2 and 102 = 1-4. Such a constel- 



lation of parameters takes the coupled system into and out of 
synchrony intermittently. The phase difference presented on 
Fig-S(a) is bounded during the synchronized intervals lEoll . 
Applying the procedure of return maps and the modified New- 
ton root-finding method [Q] we detected the synchronization 
intervals (Fig. [4] (b)) consistently with the phase difference. 
The map procedure is equivalent to determining if a model 
of coupled phase oscillators, with parameters as the one in- 
ferred, is synchronized. Note also that during the synchro- 
nized intervals, the phases become almost identical and they 
do not span enough the available space for inference. This can 
result in inferred parameters that are far from the intrinsic val- 
ues. However, the whole set of parameters is again correlated 
as if it was coming from a synchronized system. One might 
ask - why we need to construct such a 'complicated' proce- 
dure to detect synchronization, when something as simple as 
phase difference can give the similar answer. The point is that 
with the use of the intrinsic inferred parameters one can distin- 
guish whether the phase slips and synchronization transitions 
are noise-induced or not [Q]]. 

Coupling functions are arguably the most important part 
of the description of interactions. They can describe the 
functional relationship, and the law governing the mutually 
interactions and the routes to qualitative transitions. Cou- 
pling functions have been used to describe different aspects 
of the interactions of oscillatory systems of various natures, 
including cardio-respiratory, electrochemical and mechanical 



400 



'200 



(a) 



(b) 





200 



400 600 
Time [s] 



800 



1000 



Figure 4: Detection of intermittent synchronization from system l[8j. 
Phase difference ip(t) showing the qualitative statistical occurrence 
of synchronized intervals - these appeared as bounded plateaus (a). 
Synchronization index demonstrating the detected intrinsic synchro- 
nization intervals (b), evaluated from the inferred parameters. High 
values J sync = 1 denote the synchronized intervals. 
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Figure 5: Coupling functions from the model of two interacting limit-cycle oscillators ([8}. The functional influence qi (0i, fa) from the second 
to the first oscillator is shown in (a); similarly the coupling function qi(<f>i, <f>2) from the first to the second oscillator is shown in (b). The 
functions are inferred within one window of data. 



By representing the phase dynamics on a 2-k- 
phase grid evaluated for the relevant inferred parameters, we 
can effectively determine and visualize them. The coupling 
functions of systems ((8]) for one window are shown in Fig. [5] 
By inferring the dynamics while allowing the time-variability 
to be traced, we can follow the actual time-evolution of the 
functions 01. Moreover, the coupling function qi((f>i, 4>j) can 
be further decomposed on self, direct and indirect coupling 
influences i23ll . and each one can be studied separately. 

V. GENERALIZATIONS 

The aim of the discussion in the previous sections was to 
present the implementation of the Bayesian inference in a 
clear and simple way. Needless to say, the framework is much 
broader and a number of important generalizations are possi- 
ble. 

The presentation mainly concentrated on the inference of 
phase dynamics. The main reason for this was to present a 
method that will be generally applicable to coupled oscilla- 
tory systems. However, there are many situations where the 
dynamics need to be analyzed directly from the measured sig- 
nals in state space. For example, the estimation of the phases 
from chaotic systems can be problematic, while the inference 
in the state domain is not at all. If state signals are about to be 
analyzed, then Eq. (Q} will have different type of base func- 
tions, e.g. polynomial, while all the rest of the equations and 
algorithms can be equally applied on the state signals. For ex- 
amples of inference of (coupled) oscillatory dynamics in state 
space see |3|-|5|]. 

The examples here included only two coupled systems, 
while in general the technique can be applied to larger group 
i.e. small-scale network of oscillators H. The phase decom- 
positions can be applied for pairwise couplings but, more im- 
portantly, joint coupling influences can also be inferred. In 



such cases, the effective coupling lf24ll can be distinguished 
from the structural by appropriate use of surrogate testing 
JH. 

The procedure of information propagation that allows time- 
variability to be tracked, depends on p w which acts as a free 
parameter. One can further improve this by making the pa- 
rameter adaptive in order to follow the time-variability more 
closely and to infer the noise more precisely. This might be re- 
alized by determining the optimal parameter from a spectrum 
of values within each window. Even though this procedure 
will be very slow, it might prove helpful and needed in some 
cases. 

Within the framework of the phase dynamics inference, the 
phase differential (j>„ is calculated as <fr n = (4>„+i — 4>n)/h. 
Improved performance can be accomplished if one provides 
this as directly estimated instantaneous frequency - e.g. us- 
ing synchrosqueezed |f2lj or nonlinear mode decomposition 
112611 . The phase base functions are not strictly restricted to 
be from the Fourier order decomposition, and other additional 
functions can be included. For example, if we have expan- 
sion up to second order K = 2, one can include also other 
components such as sin(7</>i — 4^2) in order to detect syn- 
chronization more precisely for 7:4 synchronization ratios. 



VI. CONCLUDING REMARKS 

We presented a tutorial intended to familiarize the reader 
with a technique for Bayesian inference of time-evolving cou- 
pled dynamics in the presence of noise. The tutorial gives a 
relatively comprehensive description of the method, includ- 
ing theoretical constraints, algorithms, implementation and 
demonstration of the main features on few characteristic ex- 
amples. MatLab codes for the method and the examples are 
also available. This tutorial could lead the readers to new in- 
sights and little researched phenomena, promising to be a use- 
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ful tool for signal processing of a diverse nature. 
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